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We solve the (2 + 1)D nonlinear Helmholtz equation (NLH) for input beams that collapse in the 
simpler NLS model. Thereby, we provide the first ever numerical evidence that nonparaxiality and 
backscattering can arrest the collapse. We also solve the (1 + 1)D NLH and show that solitons with 
radius of only half the wavelength can propagate over forty diffraction lengths with no distortions. 
In both cases we calculate the backscattered field, which has not been done previously. Finally, 
we compute the dynamics of counter-propagating solitons using the NLH model, which is more 
comprehensive than the previously used coupled NLS model. 

The nonlinear Schrodinger equation (NLS) is the canonical model in nonlinear optics for propagation of intense 
laser beams in isotropic Kerr media. In the case of propagation through a bulk medium, Kelley [l| used the 2D 
NLS to predict the possibility of a catastrophic collapse of beams whose input power is above the critical power for 
collapse. In the case of propagation through planar waveguides, the ID NLS was used to predict the existence of 
spatial solitons H. Both beam collapse in bulk medium and spatial solitons in planar waveguides were observed in 
experiments d, 0]. More recently, configurations of two counter-propagating beams were modeled by two coupled 
NLS equations [5[. 

In nonlinear optics, the NLS is derived from the nonlinear Maxwell equations via a series of approximations. First, 
if the electric field is monochromatic and third harmonic generation is neglected, Maxwell's equations reduce to the 
vector nonlinear Helmholtz equation (NLH). If the field is also linearly polarized, the vector NLH reduces to the scalar 
NLH [6] . Finally, the NLS is derived from the scalar NLH using the paraxial approximation, which is valid when the 
beam radius is sufficiently large compared with the wavelength. As, however, the 2D NLS predicts that the beam 
radius shrinks to zero at collapse, the paraxial approximation breaks down at this point. In the case of spatial ID 
solitons, the paraxial approximation sets a lower limit on the soliton radius. 

The singular behavior of the 2D NLS solutions for collapsing beams is non-physical. Therefore, an important 
question is whether the singularity formation is already arrested by taking one step back in the aforementioned 
series of approximations and employing the scalar NLH model, or only in a more comprehensive model. Both the 
mathematical analysis and simulations of the scalar NLH have proved to be considerably more difficult than for 
the NLS, since for the NLH one solves a nonlinear boundary- value problem, whereas the NLS requires solving an 
initial value problem. An additional computational obstacle is that unlike the NLS, which governs the slowly varying 
envelope, the NLH has to be approximated with sub- wavelength resolution. For these reasons, the question of collapse 
in the scalar NLH model was not fully answered for over 40 years. 

Previously, numerical simulations and asymptotic analysis [3, [H, Q suggested that nonparaxiality arrests the collapse 
in bulk medium. These studies, however, applied various simplifying approximations to the scalar NLH. In partic- 
ular, they considered only forward traveling waves and completely neglected the backscattered field. Even though 
backscattering is generally believed to be "small", it may still significant affect the overall propagation, because 
collapse dynamics in the 2D cubic NLS is extremely sensitive to small perturbations [Io| . 

To study the arrest of collapse in the scalar NLH with no simplifying assumptions (and in particular, with the 
backscattering included), Fibich and Tsynkov developed a fixed-point iterative numerical method for solving the NLH 
as genuine boundary value problem [ll|,ll2|, which is based on freezing the nonlinearity at each iteration. This method 
converged for input powers below the critical power for collapse P cr , but diverged for input powers higher than P cr . It 
was unclear, however, whether the divergence above P cr was due to limitations of the numerical method, or because 
collapse is not arrested in the scalar NLH model. Subsequently, the method of [ll|, [l2[ was used to show numerically 
the arrest of collapse in the linearly-damped scalar NLH [13]. In that work, however, the magnitude of damping was 
much larger than in actual physical settings, and could not be reduced to zero. More recently, Sever proved existence 
of solutions (and hence arrest of collapse) in the scalar NLH with self-adjoint boundary conditions [141 ]. The proof 
in [Hj, however, relies heavily on self- adjoint ness, whereas propagating fields satisfy radiation boundary conditions 
(BCs), which are non self-adjoint. Therefore, until now, there has been no conclusive evidence that the collapse is 
arrested in the scalar NLH model. 

In [TEj . we studied numerically the (0 + 1)D NLH, which models the propagation of plane waves in a Kerr medium. 
In this case, the solution always exists, but becomes non-unique (bistable) above a certain input power threshold [lr| . 
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FIG. 1: The physical setup: A: single beam, B: counter-propagating beams. C: A schematic of the upstream BC at z — Z max + £, 
which freely admits all forward propagating waves (red). D: A schematic of the downstream BC at z — —8, which freely admits 
all the backward propagating waves (blue) to pass, and also specifies the (forward moving) incoming beam. 



Numerically, we observed that the fixed-point frozen nonlinearity method of pj],[l2[ converges for low input powers, but 
diverges for higher powers which are still below the threshold for non- uniqueness. This indicates that the divergence of 
the fixed-point frozen nonlinearity method is due to the numerical methodology itself, rather than to non-uniqueness or 
non-existence of the solutions. Therefore, an alternative iterative solver, based on Newton's method, was constructed 
and shown to have much better convergence properties. In this Letter, we extend the Newton-based method of [HI 
to the multi-dimensional case. The resulting technique enables us to solve the (2 + 1)D NLH for input powers above 
P cr . Hence, we obtain the first ever computational evidence that the collapse of the beam is indeed arrested in the 
scalar NLH model. We also calculate the field backscattered from the domain. Moreover, we solve the (1 + 1)D NLH 
for a "nonparaxial" soliton with radius equal to half a wavelength, and observe that it propagates virtually unchanged 
over 40 diffraction lengths. This indicates that such beams are still in the paraxial regime. Finally, we solve the 
(1 + 1)D NLH for two counter-propagating beams and compare the results to those obtained using the coupled NLS 
model. 

The propagation of linearly polarized, continuous wave beams in isotropic Kerr media is governed by the scalar 
nonlinear Helmholtz equation: 

E zz (z, x ± ) + A ± E + k 2 (1 + (2n 2 /n )\E\ 2 ) E = 0, (1) 

where E is the electric field, fco is the linear wavenumber, no is the linear index of refraction and n 2 is the Kerr 
coefficient. In the bulk medium (2 + 1)D case xx = (x,y) and Ax = d 2 + d 2 ] in the planar waveguide (1 + 1)D 
case xx = x and Ax = d 2 . We consider an incoming beam traveling in the positive z direction (henceforth "forward" 
or "right") impinging on a finite- length Kerr material slab at the z = interface and exiting the Kerr medium 
at the z = Z max interface, see Fig. DJA). A portion of the field may be reflected by the interfaces at z = or 
z = ^max, or backscattered inside the Kerr medium, because of the variations of the index of refraction induced by 
the forward-propagating beam. To derive the NLS, the standard approach is to represent the field as E = Ae lk ° z , 
where the envelope A is assumed slowly varying. Using the standard rescaling xx = xx/^o, z = z/2Ldf and 
A(z,5t±) = 1/2^2/^0^0^0 • ^(2, xx), where ro is the input beam radius and Ldf = ^0^0 * s ^ ne diffraction length, the 
NLH can be written in the dimensionless form 

f 2 A zz (z, xx) + iA~ z + Axi + \A\ 2 A = 0, (2) 

where f 2 — (r§k§)~ 2 = (2^) 2 * s the nonparaxiality parameter. Typically Ao <C ro so that f 2 <C 1 and f 2 A zz <C A z . 
Therefore, the paraxial approximation, which consists of neglecting f 2 A zz , leads to the NLS 

iA 2 (z, xx) + Axi + \A\ 2 A = 0. (3) 

In our simulations, the (2 + 1)D NLH with cylindrical symmetry, i.e., E = E(z,r) where r = |xx| = y 1 x 2 + y 2 , 
is approximated with a fourth order finite- difference scheme. The solution is computed for — 5 < z < Z max + 6 in 
order to implement the BCs in the linear regions. At the material interfaces z = and z = Z max where the index of 
refraction is discontinuous, Maxwell equations for a normal-incident field imply that E and E z are continuous across 
the interfaces. At z = Z max + S we imposed the radiation BC that the field does not have any left-going component 
for z > Z max , see Fig. [Tp. Similarly, at z = —5 we implement the two-way radiation BC that for z < the field 
does not have right-going components except for the prescribed incoming beam which impinges on the interface z = 
with a transverse profile E{ nc (r), see Fig. [Tp. Because z = — 5 and z = Z ma x + S are outside the Kerr slab, the field 
propagation there is linear, which simplifies the implementation of the radiation BCs, se e ITU [l2| for more details. 
The discretized system of nonlinear algebraic equations is solved using Newton's method jl5j ]. 

In order to focus on the effects of the Kerr nonlinearity, the values of no in the Kerr medium (0 < z < Z max ) and 
in the surrounding linear medium (z < and at z > Z max ) are chosen to be equal, so that to eliminate the reflections 
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FIG. 2: (color online) Arrest of collapse in the (2 + l)D NLH. A: \E\ 2 . B: S z . 



due to discontinuity of no at the interfaces. However, discontinuities in the nonlinear coefficient are not eliminated, 
and are a source of reflections at z = 0, Z max . Our numerical method can be applied to the case of different no with 
no change [15]. Note that, since we solve the NLH in non-dimensional form, simulations in this Letter are valid for 
any physical value of ko , no , ri2 that corresponds to the same dimensionless quantities f 2 and P/P cr . 

-l 



The (2 + 1)D NLH (pQ) was solved for an incoming collimated Gaussian beam E- inc = ( y^W^o^o^o) J e 



-(r/ro) 2 

of radius ro = 1.27Ao, corresponding to nonparaxiality parameter of f 2 — (fco^o) -2 = 1/64, and input power of 
P = 1.29P cr . The NLH solution initially self-focuses, until z « 0.8Ljjf where the collapse is arrested, after which 
the solution defocuses, see Fig. EfA). The corresponding NLS solution collapses at z c = 0.68L^i?, see Fig. [3j This 
comparison of the NLH and NLS provides a direct numerical evidence that collapse is arrested in the scalar NLH 
model. 

The fast oscillations of \E\ 2 in the z direction in Fig. [2](A) are 
not a numerical artifact, but rather account for the actual physics. 
Indeed, let us first note that a part of the forward-propagating 
wave is reflected backwards by the material interfaces at z = 
and z = Z max . In addition, since the forward propagating beam 
induces changes in the refraction index, part of the beam may 
be backscattered inside the Kerr medium. The presence of both 
forward and backward traveling fields, i.e, 



Ae 



ikoz 



Be 



ik z 



(4) 



implies that \E\ 2 « \A\ 2 + \B\ 2 + 2Re (AB*e i2/e ° 2 ) . Hence, \E\ 2 
should undergo oscillations with wave number ~ 2fco- Note that 
the analytical solutions of the (0 + 1)D NLH also exhibit these 
2/cq intensity oscillations [16[. The prediction that the intensity 
undergoes 2ko oscillations implies that the index of refraction also 
oscillates. In other words, the backward traveling field induces a 




FIG. 3: Arrest of collapse in the (2+ 1)D NLH: 
comparison of normalized on- axis \E\ 2 (blue solid), 
S z (red dashed), and NLS solution (black dotted). 



4 




FIG. 4: (color online) NLH solutions with ro/Ao = § (blue, dots), ^ (red, dash) and |- (green, solid). Solid black line is the NLS 
solution. A: Normalized on-axis intensity \E/E(z = 0)| 2 . B: Normalized on-axis Poynting flux S z /S z (z = 0). C: Transverse 
profile of the backward field at z = 0— . 




FIG. 5: (1+1)D NLH soliton with r = 4f propagating over 40L DF . A: S z . B: On-axis |£| 2 . 



2ko Bragg grating. This prediction may be tested by pump-probe 
experiments. 

In order to find a smoother representation of the solution, recall that for the NLS (|3|) the conserved beam power 
is Pnls = J\A\ 2 d5t±. For the NLH (pQ), however, the conserved beam power is Pnlh = J S z dx±, where S = 
ko!m(E*\7E) is the energy flux, or Poynting vector, and S z = fcolm (E*^-) is its z-component. Specifically, for the 
field (|4|) the value of S z reduces to the flux difference S z « fc§ (\A\ 2 — \B\ 2 ). It is therefore much smoother than \E\ 2 , 
and provides a "more natural" depiction of the NLH solution, as confirmed by comparing S z of Fig. EfB) with \E\ 2 
of Fig. [2 A). The energy flux S z shows the arrest of collapse and the focusing-defocusing dynamics more clearly, see 
also Fig. 03 

In order to analyze the effect of the nonparaxiality parameter / 2 , in Fig. [H we fix the wavelength and vary the input 
beam radius ro (while keeping the power unchanged) so that f~ 2 = 36, 64, and 144. All the NLH solutions initially 
follow the collapsing NLS solution, but later the collapse is arrested and the solution defocuses. As expected, for a 
wider input beam (lower nonparaxiality), the deviations from the NLS solution and the arrest of collapse occur later, 
and the maximum self- focusing is higher. Again we see that \E\ 2 has 2ko oscillations (whose magnitude increases as 
the input beam becomes more nonparaxial) , while the energy flux S z is smooth. 

Our numerical algorithm for solving the NLH also enables the computation of backscattering from the Kerr medium. 
In Fig.[H(C), we present the backward propagating field profiles for the previous three solutions, just before the material 
interface at z = 0. To the best of our knowledge, this is the first ever calculation of the backs cattered field of collapsing 
beams, which is due to backscattering from inside the Kerr medium and reflections from the nonlinear interface. As 
the input beam radius ro decreases, the power of the backscattered field increases from 0.46% to 0.63%, to 2.1% of 
the incoming beam power. This, as well as a comparison of magnitudes of oscillations in Fig. [H^A), shows that the 
backscattered field increases as the input beam becomes more nonparaxial. 

In (1 + 1)D configurations, the NLS possesses stable soliton solutions. It is generally believed that the paraxial 
approximation breaks down when the beam width becomes comparable to Ao and that, therefore, no solitons of 
such narrow width exist. To see that this is not the case, the (1 + 1)D NLH (pQ) is solved for the incoming NLS- 
soliton profile E- inc = (^2n2 / noro&o) -1 sech(x/ro) with width ro = Ao/2, impinging on a Kerr slab of finite length 
Z max = AOLjjf. As in the (2 + 1)D case, we impose continuity of E and E z at the material interfaces z = and 
z — ^max, and apply the radiation BCs in the linear regions at z = — S and z = Z mSiX -\-S. The solution inside the Kerr- 
slab resembles a "nonparaxial soliton" which propagates virtually unchanged, see Fig. [U[A). We note that even for 
such a narrow beam, the nonparaxiality parameter is still moderate, as f 2 = 1/tt 2 ~ 0.1, which may explain why there 
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FIG. 6: (color online) Energy flux (S z ) of the (1 + 1)D NLH with counter-propagating beams. A: Positive (forward) flux is 
red, negative (backward) flux is blue. B: The right-going beam at its incoming (blue dashed) and outgoing interface (red solid). 
Green dotted line is the coupled-NLS solution at the outgoing interface. 



there still exist soliton-like solutions. Similarly to the (2 + 1)D case, because a part of the forward propagating beam 
is backscattered, \E\ 2 exhibits the fast 2/co oscillations (Fig.EfB)), while S z is smooth. In this case, the backscattered 
field leads to 10% oscillations in \E\ 2 . 

Posada, McDonald and New [l7|,LL8| studied solutions of the (1 + 1)D Helmholtz equation over a semi-infinite Kerr 
medium, of the form A(z, x) = (W2ri2/noroko)~ 1 sech(x/ro) • e lc ' z . In later works they found similar stationary states 
for different nonlinear it ies pjl [20[. These solutions do not have any backward propagating components. In contrast, 
for the finite- length Kerr medium simulation of Fig [51 some backward moving waves must exist, because of reflections 
from the material discontinuity at z = Z max , and the full NLH as a boundary- value problem must be solved. 

Another (1 + 1)D configuration of recent interest is that of counter-propagating beams, when a right traveling 
soliton impinges at the left interface and a left traveling beam impinges at the right interface (Fig. [2(B)). This 
configuration was analyzed numerically by Cohen et al.[5] using a coupled NLS system, which is derived from the 
NLH by employing the paraxial approximation and further assuming that asynchronous terms of the Kerr nonlinearity 
can be neglected. In doing so, the BCs should simultaneously account for the coupled incoming and outgoing fields 
at each interface. As noted in these BCs can only be approximately accommodated in the coupled NLS model. In 
contrast, they can be fully implemented in the NLH model, without any approximation. Fig. [6] presents our solution 
of the NLH for counter-propagating beams of radius ro = Ao that enter a Kerr material slab at the opposite interfaces 
with a transverse displacement of d = 4.4Ao, and propagate over IOLdf- It shows that the beams are slightly attracted 
toward each other and also become wider as they propagate. The results are in close agreement with the coupled 
NLS model, see Fig. [61(B). Therefore, the more comprehensive NLH model confirms the validity of the coupled NLS 
model for counter-propagating beams even for the "extreme" parameters of ro = Ao and d = 4.4Ao- 

In this work, we solve the scalar NLH (pp), which is the simplest model for the propagation of light in Kerr media 
that incorporates nonpar axial effects, backscattering and reflections from material interfaces. Moreover, unlike the 
NLS, it accurately models the propagation of oblique beams and reflections from interfaces at arbitrary angles. This 
model neglects vectorial effects, i.e., linear and nonlinear coupling between the three components of the electric field. 
We note that the vectorial effects scale as / 2 , and hence are of the same order of magnitude as nonparaxiality. In 
bulk media, they have been shown to have the same effect as nonparaxiality, which is to arrest the collapse [y, [2l| . 
In contrast, in planar waveguides, the solitons are stable. Hence, when ro = Ao/2, f 2 is small, and so we expect that 
both nonpar axial and vectorial effects are likely to have a secondary effect on the propagation dynamics. Therefore, 
the sub-wavelength solitons predicted in the Letter are likely to remain stable also in the more comprehensive vector 
NLH model. 
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